function AFRIAT = HPZ_Afriat_efficiency_index (expenditure)
%UNTITLED Summary of this function goes here
%   Detailed explanation goes here
global Threshold

num_of_iterations = 30;

AFRIAT_UPPER = 1;

AFRIAT_LOWER = 0;

AFRIAT = 1/2;

[rows,cols] = size(expenditure);

for i=1:num_of_iterations
    
    af_exp = expenditure - (diag(diag(expenditure))*(1-AFRIAT));
   
    %The matrix REF has at the cell in the i'th row and the j'th
    %column, the difference between the value of the bundle that was chosen in 
    %observation i and the bundle that was chosen in observation j given the 
    %prices of observation i

    REF = diag(af_exp)*ones(rows,1)' - af_exp;

    %The matrix DRP has at the cell in the i'th row and the j'th
    %column, 1 if and only if the bundle that was chosen in 
    %observation i is directly reveal prefered to the bundle that was chosen 
    %in observation j (the corresponding value in REF is greater than -1).

    DRP = ceil((REF+Threshold)/(max(max(abs(REF+Threshold)))+1));

    %The matrix DRP has at the cell in the i'th row and the j'th
    %column, 1 if and only if the bundle that was chosen in 
    %observation i is directly reveal prefered to the bundle that was chosen 
    %in observation j (the corresponding value in REF is greater than -1).

    SDRP = ceil((REF-Threshold)/(max(max(abs(REF-Threshold)))+1));

    % statement needed for the graph theory external package

    set_matlab_bgl_default(struct('full2sparse',1));

    %The matrix NS_RP has at the cell in the i'th row and the j'th
    %column, Inf if and only if the bundle that was chosen in 
    %observation i is not reveal prefered to the bundle that was chosen 
    %in observation j. Otherwise it includes a positive integer.

    NS_RP = all_shortest_paths(DRP);

    %Create RP

    RP = eye(rows,cols);

    %The matrix RP has at the cell in the i'th row and the j'th
    %column, 1 if and only if the bundle that was chosen in 
    %observation i is reveal prefered to the bundle that was chosen 
    %in observation j. 

    for j=1:rows
        for k=1:cols
            if ~isinf(NS_RP(j,k)) && ~(k==j)
                RP(j,k)=1;
            end
        end
    end

    %To test for GARP we will do the following: for every pair of
    %choices x and y if xRy then not yP0x. We will take RP and the transpose of
    %SDRP and multiply element by element. Every 1 correspondes to 
    %xRy and yP0x. The final matrix is the zero matrix if and only if 
    %GARP is satisfied. 

    GARP = RP.*(SDRP');

    GARP_ERRORS = sum(sum(GARP));

    %If the GARP matrix is only zeros then the data satisfies GARP
    if GARP_ERRORS==0
        AFRIAT_LOWER = AFRIAT;
    else
        AFRIAT_UPPER = AFRIAT;
    end 

    AFRIAT = (1/2)*(AFRIAT_LOWER + AFRIAT_UPPER);        
end


















end

